Hands-On Exercise 03

Author

Henry Low

Published

September 7, 2024

Modified

September 18, 2024

Data Sources

(saved under ‘data’ folder)
Study Area - Punggol Planning Area
Punggol_St - line feature geospatial dataset containing road network
Punggol_CC - point feature geospatial dataset containing location of childcare centres

Chapter 7: Network Constrained Spatial Point Patterns Analysis

7.1 Setting Up

7.1.1 Loading the R packages

pacman::p_load(sf, spNetwork, tmap, tidyverse)

7.1.2 Importing spatial data

# Load network dataset
network <- st_read(dsn="data", layer="Punggol_St") # Check that geometry is Linestring and not multilinestring
Reading layer `Punggol_St' from data source 
  `C:\Users\Henry\Desktop\SMU Masters\2024-2025 T1\Geospatial Analytics & Applications\Project\GeospatialWebsite\Hands-On_Ex\Hands-On_Ex_03\data' 
  using driver `ESRI Shapefile'
Simple feature collection with 2642 features and 2 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 34038.56 ymin: 40941.11 xmax: 38882.85 ymax: 44801.27
Projected CRS: SVY21 / Singapore TM
# Load childcare dataset
childcare <- st_read(dsn="data", layer="Punggol_CC") %>%
  st_zm(drop = TRUE, what = "ZM") # Remove Z dimension
Reading layer `Punggol_CC' from data source 
  `C:\Users\Henry\Desktop\SMU Masters\2024-2025 T1\Geospatial Analytics & Applications\Project\GeospatialWebsite\Hands-On_Ex\Hands-On_Ex_03\data' 
  using driver `ESRI Shapefile'
Simple feature collection with 61 features and 1 field
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 34423.98 ymin: 41503.6 xmax: 37619.47 ymax: 44685.77
z_range:       zmin: 0 zmax: 0
Projected CRS: SVY21 / Singapore TM
# See childcare features
childcare
Simple feature collection with 61 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 34423.98 ymin: 41503.6 xmax: 37619.47 ymax: 44685.77
Projected CRS: SVY21 / Singapore TM
First 10 features:
      Name                  geometry
1   kml_10 POINT (36173.81 42550.33)
2   kml_99 POINT (36479.56 42405.21)
3  kml_100 POINT (36618.72 41989.13)
4  kml_101 POINT (36285.37 42261.42)
5  kml_122  POINT (35414.54 42625.1)
6  kml_161 POINT (36545.16 42580.09)
7  kml_172 POINT (35289.44 44083.57)
8  kml_188 POINT (36520.56 42844.74)
9  kml_205  POINT (36924.01 41503.6)
10 kml_222 POINT (37141.76 42326.36)
# See childcare crs information
st_crs(childcare)
Coordinate Reference System:
  User input: SVY21 / Singapore TM 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]
# See network features
network
Simple feature collection with 2642 features and 2 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 34038.56 ymin: 40941.11 xmax: 38882.85 ymax: 44801.27
Projected CRS: SVY21 / Singapore TM
First 10 features:
     LINK_ID                   ST_NAME                       geometry
1  116130894                PUNGGOL RD LINESTRING (36546.89 44574....
2  116130897 PONGGOL TWENTY-FOURTH AVE LINESTRING (36546.89 44574....
3  116130901   PONGGOL SEVENTEENTH AVE LINESTRING (36012.73 44154....
4  116130902   PONGGOL SEVENTEENTH AVE LINESTRING (36062.81 44197....
5  116130907           PUNGGOL CENTRAL LINESTRING (36131.85 42755....
6  116130908                PUNGGOL RD LINESTRING (36112.93 42752....
7  116130909           PUNGGOL CENTRAL LINESTRING (36127.4 42744.5...
8  116130910               PUNGGOL FLD LINESTRING (35994.98 42428....
9  116130911               PUNGGOL FLD LINESTRING (35984.97 42407....
10 116130912            EDGEFIELD PLNS LINESTRING (36200.87 42219....
# See network crs information
st_crs(network)
Coordinate Reference System:
  User input: SVY21 / Singapore TM 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]

7.2 Visualize Geospatial Data

# Visualize both geospatial data
plot(st_geometry(network)) # Use st_geometry to plot only the network without corresponding attributes
plot(childcare,add=T,col='red',pch = 19)

# Interactive Visualization with tmap 
tmap_mode('view')
tmap mode set to interactive viewing
tm_shape(childcare) +
  tm_dots() +
  tm_shape(network) +
  tm_lines()
# Set tmap_mode back to plot
tmap_mode('plot')
tmap mode set to plotting

7.3 Network KDE (NKDE) ANalysis

7.3.1 Prepare Lixel Objects

# Cut spatiallines objects into lixels with specified minimal distance using lixelize_lines from spNetwork
lixels <- lixelize_lines(network, 700, mindist = 375)
# Use nearest neighbour to gauge the appropriate distance - bottom ~25%-10%

7.3.2 Generate line centre points

# create line centre points with lines_center from spNetwork
samples <- lines_center(lixels) 

tmap_mode('view')
tmap mode set to interactive viewing
tm_shape(lixels) + 
  tm_lines() + 
  tm_shape(samples) + 
  tm_dots(size =0.01)
tmap_mode("plot")
tmap mode set to plotting

7.3.3 Perform NKDE

# # Remove Z-dimension from childcare (No longer required)
# childcare_rmz <- st_zm(childcare)

# Perform NKDE
densities <- nkde(network, 
                  events = childcare,
                  w = rep(1, nrow(childcare)),
                  samples = samples,
                  kernel_name = "quartic", # Try to avoid gaussian if intensity is negative
                  bw = 300, 
                  div= "bw", 
                  method = "simple", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg = 5, 
                  sparse = TRUE,
                  verbose = FALSE)
# Densities - represents intensity
# Insert computed density values as density field
samples$density <- densities
lixels$density <- densities

# Recaling to enhance mapping
samples$density <- samples$density*1000
lixels$density <- lixels$density*1000

# Visualize interactive plot
tmap_mode('view')
tmap mode set to interactive viewing
tm_shape(lixels)+
  tm_lines(col="density")+
tm_shape(childcare)+
  tm_dots()
tmap_mode('plot')
tmap mode set to plotting

7.4 Network Constrainted G- and K- Function Analysis

Ho: The observed spatial point events (i.e distribution of childcare centres) are uniformly distributed over a street network in Punggol Planning Area.

# Perform complete spatial randomness test using kfunction from spNetwork
kfun_childcare <- kfunctions(network, 
                             childcare,
                             start = 0, 
                             end = 1000, 
                             step = 50, 
                             width = 50, 
                             nsim = 50, # 51 simulations 
                             resolution = 50,
                             verbose = FALSE, 
                             conf_int = 0.05)
# Visualize k-function output
kfun_childcare$plotk

Conclusion
- There is non complete spatial randomness
- In certain distances (~250-500), childcares are regularly spaced out
- At higher distances (>600), childcares are randomly spaced out